rm(list=ls(all.names=TRUE))
set.seed(221269)


#### Loading data 
load("anes2008.RData")


#### Loading libraries
library(lavaan)
library(semTools)



#### Setting ordered variables ####
orderedvars <- c('respect', 'manners', 'obedience', 'behaved')



####################
#### Configural ####
####################

Config1group <- 
  '
## Factor loadings
authorit =~ c(NA)*respect + c(NA)*manners + c(NA)*obedience + c(NA)*behaved


## Thresholds
respect | c(NA)*t1
respect | c(NA)*t2

manners | c(NA)*t1
manners | c(NA)*t2

obedience | c(NA)*t1
obedience | c(NA)*t2

behaved | c(NA)*t1
behaved | c(NA)*t2


#Latent variable: mean
authorit ~ c(0)*1


#Latent variables: variances
authorit ~~ c(1)*authorit


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1)*respect
manners     ~*~ c(1)*manners
obedience   ~*~ c(1)*obedience
behaved     ~*~ c(1)*behaved


## Observed variables (latent variate): Intercepts
respect    ~ c(0)*1
manners    ~ c(0)*1
obedience  ~ c(0)*1
behaved    ~ c(0)*1
'


#### Whites ####
Config1fit.Wh <- lavaan::lavaan(
  data=anes2008, 
  model=Config1group, 
  group="race", group.label=c("White"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Config1fit.Wh, fit=T, standardized=T)
lavResiduals(Config1fit.Wh)
# fitMeasures(Config1fit.Wh)
modificationIndices(Config1fit.Wh)

# Observed corr
lavInspect(Config1fit.Wh, what = "obs")$cov
# Model-implied corr
lavInspect(Config1fit.Wh, what = "implied")$cov


#### Blacks ####
Config1fit.Bl <- 
  lavaan::lavaan(
    data=anes2008, 
    model=Config1group, 
    group="race", group.label=c("Black"), 
    ordered=orderedvars, 
    meanstructure=TRUE, parameterization="delta", 
    test="scaled.shifted", se="robust.sem", estimator="ULS",
    mimic="lavaan", model.type="cfa", start="simple", 
    # Weights
    sampling.weights="weight", 
    sampling.weights.normalization="group", 
    # Identification constraints
    int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
    std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
    auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
    auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
    # Categorical estimation
    zero.cell.warn=TRUE, 
    zero.add=c(0.5,0.5), 
    zero.keep.margins=TRUE, 
    # Computational / optimization options
    check.start=TRUE, check.post=TRUE,
    check.gradient=TRUE, check.vcov=TRUE, 
    check.lv.names=FALSE,
    do.fit=TRUE, optim.method="nlminb",
    implied=TRUE, verbose=TRUE)

summary(Config1fit.Bl, fit=T, standardized=T)
lavResiduals(Config1fit.Bl)
# fitMeasures(Config1fit.Bl)
modificationIndices(Config1fit.Bl)




#### Residual corr added ####

Config2group <- 
  '
## Factor loadings
authorit =~ c(NA)*respect + c(NA)*manners + c(NA)*obedience + c(NA)*behaved


## Thresholds
respect | c(NA)*t1
respect | c(NA)*t2

manners | c(NA)*t1
manners | c(NA)*t2

obedience | c(NA)*t1
obedience | c(NA)*t2

behaved | c(NA)*t1
behaved | c(NA)*t2


#Latent variable: mean
authorit ~ c(0)*1


#Latent variables: variances
authorit ~~ c(1)*authorit


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1)*respect
manners     ~*~ c(1)*manners
obedience   ~*~ c(1)*obedience
behaved     ~*~ c(1)*behaved


## Observed variables (latent variate): Intercepts
respect    ~ c(0)*1
manners    ~ c(0)*1
obedience  ~ c(0)*1
behaved    ~ c(0)*1


## Residual correlation, freely est
respect ~~ c(NA)*behaved
'


#### Whites ####
Config2fit.Wh <- 
  lavaan::lavaan(
    data=anes2008, 
    model=Config2group, 
    group="race", group.label=c("White"), 
    ordered=orderedvars, 
    meanstructure=TRUE, parameterization="delta", 
    test="scaled.shifted", se="robust.sem", estimator="ULS",
    mimic="lavaan", model.type="cfa", start="simple", 
    # Weights
    sampling.weights="weight", 
    sampling.weights.normalization="group", 
    # Identification constraints
    int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
    std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
    auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
    auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
    # Categorical estimation
    zero.cell.warn=TRUE, 
    zero.add=c(0.5,0.5), 
    zero.keep.margins=TRUE, 
    # Computational / optimization options
    check.start=TRUE, check.post=TRUE,
    check.gradient=TRUE, check.vcov=TRUE, 
    check.lv.names=FALSE,
    do.fit=TRUE, optim.method="nlminb",
    implied=TRUE, verbose=TRUE)

summary(Config2fit.Wh, fit=T, standardized=T)
lavResiduals(Config2fit.Wh)
# fitMeasures(Config2fit.Wh)
modificationIndices(Config2fit.Wh)


#### Blacks ####
Config2fit.Bl <- 
  lavaan::lavaan(
    data=anes2008, 
    model=Config2group, 
    group="race", group.label=c("Black"), 
    ordered=orderedvars, 
    meanstructure=TRUE, parameterization="delta", 
    test="scaled.shifted", se="robust.sem", estimator="ULS",
    mimic="lavaan", model.type="cfa", start="simple", 
    # Weights
    sampling.weights="weight", 
    sampling.weights.normalization="group", 
    # Identification constraints
    int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
    std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
    auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
    auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
    # Categorical estimation
    zero.cell.warn=TRUE, 
    zero.add=c(0.5,0.5), 
    zero.keep.margins=TRUE, 
    # Computational / optimization options
    check.start=TRUE, check.post=TRUE,
    check.gradient=TRUE, check.vcov=TRUE, 
    check.lv.names=FALSE,
    do.fit=TRUE, optim.method="nlminb",
    implied=TRUE, verbose=TRUE)

summary(Config2fit.Bl, fit=T, standardized=T)
lavResiduals(Config2fit.Bl)
# fitMeasures(Config2fit.Bl)
modificationIndices(Config2fit.Bl)



#### Model improvement with add of resid corr

## Whites
lavTestLRT(Config1fit.Wh, Config2fit.Wh, method="satorra.bentler.2001")

## Blacks
lavTestLRT(Config1fit.Bl, Config2fit.Bl, method="satorra.bentler.2001")
